#include <iostream>
#include <cmath>
#include <algorithm>
#include <stdio.h>
#include <stdlib.h>
#include <fstream>
#include <sys/time.h>
#define epsilon 1.e-8
using namespace std;
#include "mpi.h"


int main(int argc, char **argv) {

  int myrank;
  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);

  int M,N;
  int i,j,k;
  string T,P,Db;
  M = atoi(argv[2]);
  N = atoi(argv[3]);

  double elapsedTime,elapsedTime2;
  timeval start,end,end2;

  if(argc > 4){

    T = argv[4];
    if(argc > 5){
      P = argv[5];
      if(argc > 6){
        Db = argv[6];
      }
    }
  }
  // cout<<"M-->"<<M<<"   N--->"<<N<<endl;
 // // cout<<T<<P<<endl;  
  
  double **U_t;
  double alpha, beta, gamma,**Alphas,**Betas,**Gammas;
  int acum = 0;
  int temp1, temp2;
 

  U_t = new double*[M];  //二维数组的动态分配是否有问题
  Alphas = new double*[M];
  Betas = new double*[M];
  Gammas = new double*[M];

  for(int i =0; i<M; i++){
      U_t[i] = new double[N];
      Alphas[i] = new double[M];
      Betas[i] = new double[M];
      Gammas[i] = new double[M];
   }

  ifstream matrixfile("matrix");
  if(!(matrixfile.is_open())){
    cout<<"Error: file not found"<<endl;
    return 0;
  }

  for(i = 0; i < M; i++){
    for(j =0; j < N; j++){
      matrixfile >> U_t[i][j];
    }
  }
  matrixfile.close();

  int nprocs;
  MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
  gettimeofday(&start, NULL);
  int nSize = M/nprocs;
  int m_start = nSize*myrank;
  int m_end  = nSize*(myrank+1);
  if(myrank > 0){
    m_start = m_start + M%nprocs;
  }
  m_end = m_end+M%nprocs;
  // cout<<"m_start "<<m_start<<" m_end "<<m_end<<"   nSize: "<<nSize<<endl;
  
   for(i =m_start; i<m_end;i++){    //convergence

    for(j = 0; j<M; j++){

      alpha =0.0;
      beta = 0.0;
      gamma = 0.0;
      for(k = 0; k<N; k++){


         alpha = alpha + (U_t[i][k] * U_t[i][k]);
         beta = beta + (U_t[j][k] * U_t[j][k]);
         gamma = gamma + (U_t[i][k] * U_t[j][k]);
       }
      Alphas[i][j] = alpha;
      Betas[i][j] = beta;
      Gammas[i][j] = gamma;

     }
  }
  MPI_Barrier(MPI_COMM_WORLD);
  gettimeofday(&end, NULL);
  MPI_Status m_status;
  int len=M*nSize;
  int ind = 0;
  double *alphabuf = new double[len];
  double *betabuf = new double[len];
  double *gamabuf = new double[len];
  if(myrank ==0){  //接收数据
    int rev;
    for(rev=1;rev<nprocs;rev++){
      MPI_Recv(&alphabuf[0],len,MPI_DOUBLE,rev,rev,MPI_COMM_WORLD,&m_status);
      MPI_Recv(&betabuf[0],len,MPI_DOUBLE,rev,rev,MPI_COMM_WORLD,&m_status);
      MPI_Recv(&gamabuf[0],len,MPI_DOUBLE,rev,rev,MPI_COMM_WORLD,&m_status);
      // cout<<"myrank: "<<rev<<endl;
      // for(i=0;i<len;i++){
      //   cout<<betabuf[i]<<"  ";
      // }
      // cout<<endl;

      int mm_start = m_end+(rev-1)*nSize;
      int mm_end = m_end+nSize*rev;
      ind=0;
      // cout<<mm_start<<"--->"<<mm_end<<endl;
      for (i = mm_start; i < mm_end; i++)
      {
        for ( j = 0; j< M; j++)
        {
             Alphas[i][j] = alphabuf[ind];
             Betas[i][j] = betabuf[ind];
             Gammas[i][j] = gamabuf[ind];
             ind = ind+1;;
        }
      }
    }

  }else{  //发送数组
    ind=0;
    for(i = m_start; i < m_end; i++){
        for (j = 0; j< M; j++)
        {
             alphabuf[ind] = Alphas[i][j];
             betabuf[ind] = Betas[i][j];
             gamabuf[ind] = Gammas[i][j];
             ind++;
        }
      }
      // cout<<"myrank send "<<myrank<<endl;
      // for(i=0;i<len;i++){
      //   cout<<betabuf[i]<<" ";
      // }
      // cout<<endl;


      MPI_Send(&alphabuf[0],len,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
      MPI_Send(&betabuf[0],len,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
      MPI_Send(&gamabuf[0],len,MPI_DOUBLE,0,myrank,MPI_COMM_WORLD);
  }

// fix final result
    MPI_Finalize();

  //Output time and iterations
  if(myrank == 0)
  {

      if(T=="-t" || P =="-t"){
        elapsedTime = (end.tv_sec - start.tv_sec) * 1000.0;
        elapsedTime += (end.tv_usec - start.tv_usec) / 1000.0;
        cout<<"Time: "<<elapsedTime<<" ms."<<endl<<endl;
      }

      // Output the matrixes for debug
      if(T== "-p" || P == "-p"){
      cout<<"Alphas"<<endl<<endl;
      for(i =0; i<M; i++){

        for(j =0; j<M;j++){
              
          cout<<Alphas[i][j]<<"  ";
        }
        cout<<endl;
      }

      cout<<endl<<"Betas"<<endl<<endl;
      for(i =0; i<M; i++){

       for(j=0; j<M;j++){   
          cout<<Betas[i][j]<<"  ";
       }
       cout<<endl;
      }

      cout<<endl<<"Gammas"<<endl<<endl;
      for(i =0; i<M; i++){
        for(j =0; j<M; j++){

           cout<<Gammas[i][j]<<"  ";
      
         }
        cout<<endl;
      }

      }

      //Generate files for debug purpouse
       if(Db == "-d" || T == "-d" || P == "-d"){
        ofstream Af;
        //file for Matrix A
        Af.open("AlphasMPI.mat"); 


       // Af<<M<<"  "<<N;
        for(i = 0; i<M;i++){
          for(j =0; j<M;j++){
            Af<<" "<<Alphas[i][j];
          }
          Af<<"\n";
        }
        
        Af.close();

        ofstream Uf;

        //File for Matrix U
        Uf.open("BetasMPI.mat");

        
        for(i = 0; i<M;i++){
          for(j =0; j<M;j++){
            Uf<<" "<<Betas[i][j];
          }
          Uf<<"\n";
        }
        Uf.close();

        ofstream Vf;
        //File for Matrix V
        Vf.open("GammasMPI.mat");

         for(i = 0; i<M;i++){
          for(j =0; j<M;j++){
            Vf<<" "<<Gammas[i][j];
          }
          Vf<<"\n";
        }
        
        Vf.close();

        ofstream Sf;

     }
   }
   
   for(i = 0; i<M;i++){
     delete [] Alphas[i];
     delete [] U_t[i];
     delete [] Betas[i];
     delete [] Gammas[i];
     
   }
   delete [] Alphas;
   delete [] Betas;
   delete [] Gammas;
   delete [] U_t;

  return 0;


}









